Single-atom density of states of an optical lattice 
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We consider a single atom in an optical lattice, subject to a harmonic trapping potential. The 
problem is treated in the tight-binding approximation, with an extra parameter k denoting the 
strength of the harmonic trap. It is shown that the k — > limit of this problem is singular, in the 
sense that the density of states for a very shallow trap (/t — > 0) is qualitatively different from that of 
a translationally invariant lattice (k — 0). The physics of this difference is discussed, and densities 
of states and wave functions are exhibited and explained. 



^3 

o 

CD 



CO 



I 

C 

O 

o 



> 

On 

o 

m 
O 



-a 
c 

o 
o 



X 
S3 



The achievement of degeneracy in atomic gases has led 
to a surge of activity. One reason is that these systems 
provide a test-bed for theories developed in the context 
of solid-state physics. Indeed there have been experimen- 
tal observations of Bose-Einstein condensation (BEC) 1 1, 
vortex lattices 0, 01 jSuperfluid- insulator transitions 4|, 
fermion degeneracy [a| and, very recently, advances have 
been made towards Bardeen-Cooper-Schrieffer (BCS) su- 
perfluidity 0. A great advantage of these systems is 
one's unprecedented control of the experimental situa- 
tion. By varying the magnetic field strength, particle- 
particle interactions can be altered. Moreover, one can 
use lasers to create an optical lattice whose parameters 
can be changed at will Mj. The latter scenario has been 
the subject of recent theoretical attention |sL IflL flfl Hll| . 

One key difference between these systems and the usual 
models of the solid state is that they are not translation- 
ally invariant, due to the trapping potential necessary to 
confine the gas. However, in the usual set-up, when the 
trap is the only external potential, the local density ap- 
proximation can be employed, and thermodynamic quan- 
tities can be simply related to their values in the absence 
of the trap 0|. Trapped atoms are, therefore, a good 
model of the continuum many-particle system — analo- 
gous, in the fermion case, to 'jellium', the ideal metal of 
many-electron theory. One thus expects that, when the 
periodic optical lattice potential is added, an experimen- 
tal model of electrons in real crystals can be obtained. 

In this letter we take a first step towards testing this 
expectation by obtaining the single-particle density of 
states (DOS). This is the starting point for understand- 
ing many key features of degenerate quantum systems. 
For example, in an ideal Bose gas it determines whether 
there is BEC > while in a weakly attractive Fermi gas 
jives the critical temperature for BCS superfluidity 
The DOS becomes a smooth curve only in the limit 
in which the trap is very weak, so that the spectrum be- 
comes a continuum. As we shall see, this limit is singular, 
in the sense that the DOS is qualitatively different from 
that of an infinite, translationally invariant lattice. 

Model. We consider a single particle, with no internal 
degrees of freedom, moving in an infinite tight-binding 
lattice with one orbital per site and hopping between 
nearest neighbours. The particle is confined to a finite 



region by a harmonic trapping potential. For simplicity, 
we assume a d-dimensional hypercubic lattice. For d > 1, 
the problem is separable, so we shall focus our attention 
on the ID case, indicating how the results for d ^ 2 follow 
as the need arises. The Hamiltonian for d = 1 is 

H = -tJ2(\j) U + M + H-c) + ~/^(aj) 2 \J) U\ , 

3 3 

(1) 

where j is a site label, a the lattice constant, k the 
strength of the trap and t the nearest-neighbour hopping 
integral. (This model also occurs in the related subject of 
the dynamical diffraction of atoms by static light fields; 
see, for example, ref. 0] and references therein ) Note 
that the trap is centred at a lattice site; see refs. [lfiuTTl] 
for a discussion of the non-trivial effect of incommensu- 
ration. The model can be generalised also by considering 
anisotropic lattices and trapping potentials. 

Before describing the properties of the model, let 
us discuss its validity. In experiments @], the lattice 
is generated by counter-propagating laser fields, giv- 
ing rise to a static, periodic potential of the form 
U(x) = -U [cos 22£ + l] , Eq. QJ is valid when the wells 
of this potential are sufficiently deep: Uq ^S> ti 2 /ma 2 , 
where m is the mass of one atom. In this intense-laser 
regime, one can approximate the orbitals at the bottom 
of each well b y th ose of a harmonic oscillator, and the 
usual method yields the tight-binding parameter as 

t = [/ exp -fyVo^/Wr 3 

smaller than the first excitation energy of the orbital. 
The effect of the overall trapping potential on t can be 
neglected provided that the potential energy difference 
between the minima of neighbouring wells, AEj, is much 
less than the barrier height, Uq. Now, AEj sa na 2 j for 
j ^> 1; hence our model is accurate for states whose 
wave functions do not extend beyond site j c = Uq/ko 2 . 
Hence the state's energy e must satisfy e <C na 2 j 2 , i.e. 
£ -C e m ax = Uq/ko 2 . We would therefore like to be sure 
that e max » t, which is true provided that Uq ^> na? . 

We are interested in the density of states (DOS), p(e) = 
J2 V $ ( e — e ^)) which is a smooth curve only in the limit in 
which the spectrum {e„} is a continuum. Let us discuss 
some limiting behaviours of p(e). The solution is well 
known for k = 0, when the lattice forms a large 'box' 



which is much 
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of length L. In this case one obtains Bloch waves with 
the band dispersion relation tk — — 2icosfca, where k is 
the wave number. When L a, ka runs continuously 
from — 7r to 7r. In this limit, pie) is non-zero only for 
|e| ^ 2t (this region we call the 'band'), and has square- 
root singularities at e — ±2t (the 'top' and 'bottom' of 
the band). For d — 2 the DOS extends from e = — it to 
it and is finite at these points, but exhibits a 'van Hove' 
logarithmic singularity at e = 

In our model, the lattice itself is taken to be infinite; 
the finite extent of the wave functions is determined in- 
stead by the strength, k, of the harmonic trap. The cor- 
responding length scale is I = y/t/~K. This suggests 
that, in this model, 'continuum limit' should be taken to 
mean I > a or, equivalently, no 2 <C t. We shall work in 
this limit henceforth (unless otherwise stated). 

Let us consider the low- and high-energy states of Q. 
For energies near the bottom of the band, the dispersion 
relation has the free-particle form efe « — 2t + h 2 k 2 /2m*, 
where m* = h 2 /2ta 2 . Thus we would expect the low- 
lying eigenstates of (QJ to resemble those of the usual 
continuum harmonic oscillator. Indeed, one can show 
explicitly that \(j) ) = N J2j e ~ aj2/2 *^ 1 \j) is the ground 
state to second order in o/L Its energy is e = —2t + 
huj* /2 + o ([a/Z] 2 ) , where uj* = n/m* . To first order in 
a/l, this equals the usual result for the harmonic oscil- 
lator, measured from e = — 2t. The low-lying states are 
obtained similarly. Since the DOS of the ID harmonic 
oscillator is constant and equal to 1/hw*, we conclude 
that the low-energy limit of the DOS of Q is given by 
p(e) -> 1/fkJ* as e ^ -2t for d = 1, For d = 2, the 
DOS of the harmonic oscillator is p (e) oc e, so we expect 
a vanishing density of states as e — > — it in this case. 

For high energies, the physical nature of the states 
is quite different. When e 3> t, the first term in Q 
can be neglected, and thus the eigenstates of the Hamil- 
tonian are the position eigenstates with energies 
ej = na 2 ] 2 /2. The high-energy ID DOS is then easily 
shown to be p(e » t) = ^2/na 2 e -1 / 2 . In 2D, the anal- 
ogous calculation gives p(e) — > 2tt / na 2 as e/t — > oo. 

WKB analysis. To extend our treatment of (QJ to 
all energies, we employ the Wentzel-Kramers-Brillouin 
(WKB) approximation it has already been ap- 

plied to the momentum-space version of the problem in 
ref. [2J. By contrast, our WKB-type analysis will be 
performed in real space 01 . It differs from traditional 
WKB in that h does not appear anywhere in Q ; instead 
we choose a as our small parameter. As in the ordinary 
WKB method, we write the wave function in the form 

^p(x) ~ exp J k(x')dx'^j . The energy is determined by 

the usual quantisation condition, § k(x) dx = 2ir (n + 7), 
where n is an integer and 7 a constant. Differentiating 
this, we obtain the DOS: p(e) = dn/de = ^rjj § k{x) dx 
(note that 7 is not required). 



On the other hand, the orbit equation is modified by 
repla cing the free-particle dispersion by the tight-binding 
form |22|: — 2tcos(ka) + ^kx 2 — e. The local wavenum- 

ber becomes k(x) = ±^ arccos ^ Krr 4 7 2e ) ■ This makes a 
significant difference, because it introduces two new turn- 
ing points. The classical turning points (where k = 0) 
still exist, and are found at x = ±ar c , where 

now, however, two new points appear for energies e > 2t. 
They are associated with Bragg reflection, and they occur 
at x = ±x fc , where 

/2e-4i „, r~e 
« fc = v f_ = ^--l. (3) 

We shall refer to them as 'Bragg turning points', and the 
region between them as 'Bragg-forbidden'. 

The quantisation condition now becomes more sub- 
tle. To see why, note that in the Bragg-forbidden region 

k(x)a — 7r + i arccosh (^ 2e ^ x ) : the wave function con- 
tinues to oscillate for \x\ < Xb, despite the damping of 
its amplitude. The quantisation condition must therefore 
include these oscillations; it becomes: 

2 J k(x) dx + = 7T (n + 7) . (4) 

Now we can obtain the DOS by differentiating J1J, re- 
membering that Xb and x c depend on e according to 12I3|I : 
we obtain a lengthy analytic expression for p(e). It is 
plotted in Fig. Efa), along with the DOS obtained from 
numerical diagonalisation of , with appropriate provi- 
sion for finite-size scaling (see below): clearly the agree- 
ment is excellent. This was foreseeable, since aside from 
a vanishing number of low-lying states, all the wave func- 
tions have slowly varying k(x) (or, in the Bragg-forbidden 
region, not varying at all); thus the WKB method is ex- 
pected to be accurate at all energies. Note that p(e) 
is qualitatively different from the DOS of the K = 
tight-binding model, while it agrees quantitatively with 
p — > l/Hu* at e = -2t and with p — * ^2/ no 2 e" 1/2 for 
e/t^> 1, the limiting expressions obtained above. 

Calculating the density of states for d > 1 is simple, 
because the Hamiltonian is separable. Consider, for ex- 
ample, the case of d = 2. It is not hard to show that 
P2d(c) = Jl^ p(e - x)p(x) d X, where p 2 D is the DOS of 
the 2D system, and p is the DOS of the ID Hamiltonian 
QJ. The result of applying this formula is presented in 
Fig. Efb), and compared with numerical results; again, 
the agreement is excellent. Note that the steep feature 
at e = and the kink at e = it both result from the log- 
arithmic singularity in the ID DOS. As e — > —it, we re- 
cover P2D 0, as expected from the analytic arguments 
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Figure 1: Numerically determined wave functions for a trap 
with Ka 2 = i/100. Each wave function has been offset along 
the y-axis by its energy, in units of the hopping integral. The 
numerical calculation uses a lattice with N = 130 sites. For 
the energies considered, there are no noticeable numerical 
artefacts. The lower and upper solid curves correspond to 
the classical and Bragg turning points, eqs. 1213 
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Figure 2: Finite-size scaling of the binned DOS for ID lat- 
tices with finite strength of the trapping potential: Ka 2 /t = 
5 ■ 1(T 6 (+), 1(T 5 (x), 1(T 4 (*) and 1(T 3 (□). In all cases 
the numerical cutoff energy, ejv, lies outside the energy 
range showed in the plot. This is illustrated in the in- 
set: it shows the DOS for fixed Ka 2 = t/1000, but calcu- 
lated using lattices with different numbers of sites: N = 
300 (□), 350 (*), 420 (x) and 2000 (+). 



above. Furthermore, the large-e asymptote is exactly the 
constant predicted by that analysis. Higher-dimensional 
densities of states may be generated in a similar way. 

Finite-size scaling: numerics. In the foregoing we ob- 
tained the DOS in the continuum limit l/a — > oo. How- 
ever in experiments l/a is finite. We should therefore 
check that our results resemble the behaviour of the sys- 
tem when it has a steeper trap. Thus we have studied the 
single-particle spectrum and wave functions by numerical 
diagonalisation of (QJ p.5]. By considering increasingly 
shallow traps, the numerical calculations also allow us to 
confirm the validity of WKB for the continuum limit. We 
have applied the same method in d — 2, starting with the 
2D version of (QJ. The results are analogous so we shall, 
as above, describe in detail the ID case only. 

In order to have a finite-dimensional Hamiltonian ma- 
trix, we must describe the lattice by a finite number of 
sites, N. This leads to numerical artefacts for e > e^v, 
defined by N = 2x c (eN)/a, as the particle visits the arti- 
ficial limits of the lattice. On the other hand we expect 
this to be essentially equivalent to the original problem 
for 6 6 7v • Calculations for increasingly large N confirm 
this (see the inset of Fig. EJ) - 

Fig. shows some numerically determined wave func- 
tions, for fixed trap strength. In spite of the relatively 
small system size (l/a = 10), the states seem well de- 
scribed by the continuum-limit solution. At low energies, 
e < 0, the wave functions resemble those of a continuum 
ID harmonic oscillator. At higher energies, the period 
of these oscillations starts to become similar to the lat- 
tice constant, and commensuration effects emerge. At 
e = 2i, as that period reaches its minimum, the Bragg 
reflection points appear. Immediately above the top of 



the band, e > 2t, the modulation of the wave number is 
evident, reaching its minimum and maximum near the 
classical and Bragg turning points, respectively. The 
rapid oscillations continue in the Bragg-forbidden region, 
where the amplitude decays exponentially. As the en- 
ergy rises further, the distance between the Bragg and 
classical turning points shrinks, so that fewer and fewer 
oscillations take place between these two points. Even- 
tually the particle is forced to localise on single sites (not 
shown) . Note that the Bragg and classical turning points 
are very accurately described by Eqs. 112131 . even for the 
relatively small value of l/a considered here. Further- 
more, the number of high-energy (localised) solutions is 
reproduced correctly. This is because the condition that 
the wave function's phase be single-valued reduces, at 
k = 7T, to the condition that the localised states be sepa- 
rated by an integer number of lattice spacings (x = na) . 

We now turn to the DOS. Obviously for finite l/a the 
spectrum is not a continuum so the DOS is not a smooth 
curve. However, one can smooth it by the following pro- 
cedure: we divide an energy range into intervals, of width 
Ae, and count the states within each interval. We then 
re-scale the vertical axis of the resulting histogram so 
that each column represents the number of states per 
unit energy, on average, in the corresponding interval. 
The result is a 'binned' DOS that looks smooth only if 
the intervals are sufficiently wide. Fig.[3shows the result 
for different values of na 2 /t. Note that we have further 
re-scaled the vertical axis by an overall factor oia/l. This 
leads to the collapse of all the data onto a single curve, 
suggesting that the continuum-limit DOS describes the 
overall distribution of energy levels even for quite steep 
traps, with only l/a ~ 100 sites. The inset illustrates 
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Figure 3: DOS of a single atom in an optical lattice for (a) 
d — 1 and (b) d = 2. Curves: continuum-limit results ob- 
tained within the WKB approximation. Histograms: numer- 
ical diagonalisation for finite-size systems with trap strengths 
K = 5 x 10~ 6 t/a 2 (a) and it = 1.5 X 1(T 2 t/a 2 (b), using 
lattices with N = 5000 and 10000 sites respectively. 

the numerical cutoff artefact mentioned above. Finally, 
Fig- El compares the DOS of two finite (but fairly large) 
systems to the WKB predictions obtained above in ID 
and 2D respectively. Evidently our WKB approach cap- 
tures the continuum limit rather well. 

Conclusion. We have calculated the DOS for a single 
atom in an optical lattice; this should be regarded as the 
logical first step towards a detailed theory of the experi- 
mentally realised many-particle systems. Our results are 
based on WKB theory, and refer to the limit 1 > a or, 
equivalently, when the trapping potential becomes flat: 
k — ► 0. Numerical diagonalisation reveals this theory to 
be extremely accurate in that case, and moreover shows 
that, for finite-size systems, the binned DOS has the same 
overall features. Our main result is that the DOS, in this 
limit, is radically different from what is obtained for a ho- 
mogeneous lattice (i.e. for k = rather than k — > 0): the 
square-root singularities in the ID case are replaced by 
a logarithmic one, and the logarithmic van Hove singu- 
larity in 2D disappears altogether. Moreover, our theory 
provides a detailed picture of how this comes about. The 
crucial feature is the new turning points, associated with 
Bragg reflection, that appear at energies above the top of 
the conduction band. The possibility of inducing Bragg 
reflection by using a time-dependent external potential 



was considered in [llj. We have shown that Bragg reflec- 
tion is, in fact, essential to understand the equilibrium 
single-particle spectrum of the optical lattice. 
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Note added. — It has been drawn to our attention that 
some of the features discussed here have recently been 
pointed out in an exact diagonalisation study [23] . 
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